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Abstract 

Alexandrov proved that any simplicial complex homeomorphic to 
a sphere with strictly non-negative Gaussian curvature at each vertex 
can be isometrically embedded uniquely in as a convex polyhe¬ 
dron. Due to the nonconstructive nature of his proof, there have yet 
to be any algorithms, that we know of, that realizes the Alexandrov 
embedding in polynomial time. Following his proof, we developed the 
adiabatic isometric mapping (AIM) algorithm. AIM uses a guided adi¬ 
abatic pull-back procedure to produce “smooth” embeddings. Tests of 
AIM applied to two different polyhedral metrics suggests that its run 
time is sub cubic with respect to the number of vertices. Although 
Alexandrov’s theorem specifically addresses the embedding of convex 
polyhedral metrics, we tested AIM on a broader class of polyhedral 
metrics that included regions of negative Gaussian curvature. One 
test was on a surface just outside the ergosphere of a Kerr black hole. 
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1 Introduction 


The problem of embedding surfaces homeomorphic to with a metric of 
positive Gaussian curvature into was posed by Herman Weyl in 1916 [^. 
A useful review and analysis of the isometric embedding of Riemannian man¬ 
ifolds in Euclidean spaces can be found in . The hrst attempt to prove the 
existence of an embedding was given by Weyl himself. He was not able to 
complete his proof, but he did make progress in outlining a solution. Follow¬ 
ing Weyl’s approach, a proof was given by H. Lewy in 1938 though his solu¬ 
tion required the components of the metric to be inhnitely differentiable [^. 
Requirements on differentiability of the metric were reduced to continuous 
fourth derivatives by L. Nirenberg in 1953 [^. They were further reduced 
to continuous third derivatives by E. Heinz in 1962 [^. In 1941, Alexandrov 
gave an approach that relied on proving the existence of a convex polyhe¬ 
dron given any convex polyhedral metric. He showed that in the limit as the 
number of vertices in the polyhedral metric goes to inhnity, one recovers the 
metric in the continuum. This result can be found in a compendium of his 
work published in 2006 [^. 

Our interest in isometric embedding stems from its application in general 
relativity. Solutions to Einstein’s equations give a four dimensional spacetime 
manifold. One can slice these manifolds in such a way to get interesting 
two surfaces such as the event horizon or ergosphere. To visualize these 
surfaces it may be best to embed them in three dimensional flat Euclidean 
space. Another application of isometric embeddings include computing quasi¬ 
local energy (QLE) in general relativity. Many dehnitions of quasi-local 
energy, such as the Brown-York and Wang-Yau energies, require isometric 
embeddings of convex surfaces in [^[^. A detailed review of the 
development of QLE is given in [^. 

Though the mathematics of global isometric embeddings of convex two- 
surfaces homeomorphic to was studied exhaustively for over 60 years, 
its numerical implementation has proven formidable and has only been ap¬ 
proached practically since the mid-1990’s. In 2006, Bondarescu et. al in¬ 
troduced a numerical method using the components of the metric in the 
continuum [^. They expanded the embedding functions using spherical 
harmonics and minimized the coefficients by optimizing the squared differ¬ 
ence in metric components. Bondarescu et. al found that their search would 
get stuck in local minimum while minimizing in the space of coefficients. 
To reduce their residuals they increased the number of coefficients for each 
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embedding function. 

In 2011, M. Jasiulek and M. Korzynski introduced an algorithm that 
closely followed the continuity method prescribed by Weyl 11 . First they 


uniformized their surface using Ricci flow, which gave them conformal rela¬ 
tions between their original metric, the round sphere metric and all inter¬ 
mediate metrics. They then embedded the spherical metric into This 
allowed them to use the embedded sphere as a starting point for solving the 
linearized embedding equations which gave them the deformation in coordi¬ 
nates such that they move from the sphere to a sufficiently close conformally 
related surface. Because they were working in the continuum, they dealt 
with typical difficulties such as hnding suitable coordinate charts, solving 
non-linear elliptical PDF’s and computing integrals and derivatives of func¬ 
tions on the surface. They were also restricted to only embedding surfaces 
with strictly positive curvature due to the inversion of the extrinsic curvature 
tensor in their procedure. 

Working in the continuum has the added complication of requiring prior 
knowledge about the sign of the Gaussian curvature of the surface. Whether 
or not the embedding equations are elliptical or hyperbolic depends on the 
sign of the curvature. Additionally one may need to provide numerical tech¬ 
niques to deal with coordinate singularities. Difficulties in the continuum are 
reasons that one might want to work in the discrete. In 1995, H-P. Nollert 
and H. Herold developed the wire-frame method which was the hrst attempt 
at hnding a discrete numerical solution to the embedding problem 12 . Their 


method used position vectors where are the internal parameters 

of the surface. Given they made a Taylor series expansion around 

points on the surface allowing them to write distances in Euclidean space 
in terms of components of the internal metric. Equating this distance with 
the actual Euclidean distance gave them a function to minimize where the 
variables are coordinates in Because their optimization function only 
restricted edge lengths, and in general polyhedrons are not uniquely deter¬ 
mined by them, their method had no means of choosing what they called 
smooth solutions. 

In 2008, D. Kane et. al wrote a pseudopolynomial time algorithm to give 
a numerical realization of Alexandrov’s embedding theorem 13 . Instead of 


following the proof by Alexandrov, they modeled their algorithm using the 
proof by A. Bobenko and I. Izmestiev 14 . Given a convex polyhedral met¬ 
ric, this method uses the variation of three dimensional curvatures within 
the interior of the polyhedron to slowly deform their surface to one that is 
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isometric. Alexandrov’s theorem says given a convex polyhedral metric there 
exists a unique convex polyhedron in Kane et. al showed that their algo¬ 
rithm hnds an approximate convex polyhedron where no edge has a dihedral 
angle greater than tt -|- e. The time taken to reach this e-convex polyhe¬ 
dron depends on several intrinsic variables of the metric and user specihed 
tolerances. Their algorithm only works for convex polyhedral metrics. 

In this paper we present the adiabatic isometric mapping (AIM) algo¬ 
rithm which is a numerical approach for embedding polyhedral metrics. Like 
D. Kane et. al’s, our algorithm produces approximately convex smooth poly¬ 
hedrons given convex polyhedral metrics. For metrics that are not convex, 
AIM produces smooth polyhedrons. The AIM algorithm borrows techniques 
from several of the algorithms mentioned above. The hrst step of AIM uni- 
formizes the initial polyhedral metric under discrete Ricci flow in an affine 
time parameter t. The second step embeds this constant curvature surface 
near the surface of a sphere. We use this embedded surface as the start¬ 
ing point for an embedding flow back to the original metric. This is in a 
similar vein to M. Jasiulek and M. Korzyhski’s algorithmm. To step from 
one conformally related polyhedron to the next, we use Newton’s method 
to minimize an objective function that depends on the edge lengths of the 
polyhedral metric and the coordinate distances in We use the coordi¬ 
nates of the constant curvature polyhedron as the initial guess for Newton’s 
method to avoid the problem of local minima found in Bondarescu et al. 
This puts us close enough to the solution so that we quickly converge to the 
global minimum. We then use the newly embedded polyhedron as the initial 
value data for our next step such that At is small enough to remain near 
the global minimum. This is repeated until we reach the original polyhedral 
metric. Taking At small is not enough to ensure that we will not produce 
non-smooth solutions as seen by Nollert and Herold. Because of this, we 
introduce the convexity routine at each time step which guides our solution 
toward smooth embeddings. Using our condition on At and the convexity 
routine, we introduce the guided adiabatic pull-back which is unique to our 
algorithm. We found that AIM is capable of embedding with accuracy in the 
edge lengths to any desired tolerance above machine precision. 
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2 Adiabatic Isometric Mapping (AIM) 

Following the approach of Alexandrov, we begin by finding a continuous 
family of polyhedral metrics Pt [0 < t < tj]. A polyhedral metric p of a 
triangulated surface is a list of its triangles by vertices {vi^Vj^v^} together 
with an assignment of a length to each edge ^ij = ViVj. Here, po is the metric 
we wish to embed and pt^ is a metric with constant Gaussian curvature 
at each vertex. Each polyhedral metric in the family pt has Nq vertices, 
Ni = 3 Nq — 6 edges and N 2 = 2Nq — 4 flat triangular faces. The squared 
edge lengths of pt are given as, 

PIM} ( 1 ) 

where the indices a and b label the vertices of pt- It should be noted that 
throughout our algorithm the structure of the triangulation of pt always 
remains the same. The structure is dehned as a list of triangles by vertices. 
After we obtain p^, we hnd a polyhedron Pt^. which is an isometric embedding 
of ptj in For Pt^ to be isometric to pt^, the coordinates for the vertices 
ofPq, 

® ^ {^a {tf) , Va {tf) , Za (t/)} , (2) 

must satisfy each of the Ni distance relations, 

^Ib (tf) = (tf) - Xb itf)f + iVa (tf) - Vb itf)f + {Za (tf) - Zb (3) 

We specify freely 6 of the 3Nq coordinates in order to mod out the trans¬ 
lations and rotational degrees of freedom. This is often done by identifying 
one of the triangles Aabc of Ptf, and hxing (1) the three coordinates of vertex 
a, (2) two of the three coordinates of vertex 6, and (3) one of the remaining 
three coordinates of vertex c. The isometric embedding problem involves a 
solution to a quadratic system of 3Nq — 6 equations and unknown coordi¬ 
nates. There are many solutions that satisfy the quadratic equations; as a 
result, one hnds oneself in a “sea of solutions” that makes solving this sys¬ 
tem of equations prohibitively difficult. Not only does one have to solve the 
non-linear sparsely coupled system, they also have to hnd a solution with the 
desired extrinsic curvature. Constrained optimization problems of this kind 
are often costly to solve. To navigate through the “sea of solutions” without 
resorting to constrained optimization, AIM uses a three step procedure: (1) 
uniformization via Ricci how that provides a dimensional reduction from 3 
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Simplicial Ricci Flow 



Solve N-2 systems 
of quadratK 
equations with 3 
equations and 
unknowns 

■f 

Anneal Step A 
■¥ 

Pop out Routine 
+ 

Anneal Step B 


Figure 1: This figure illustrates the three steps of the AIM algorithm. See. 2.1 of our 
paper, “Uniformization via Ricci Flow”, begins at pq in the upper left corner and ends 
at ptf in the upper right corner. See. 2.2, “Uniform Surface Embedding” begins at ptj 
and ends at the bottom right of the figure. Finally See. 2.3, “Guided Adiabatic Pull-Back 
(GAPGAP)”, starts at Ptf in the lower right corner and ends at Pq in the lower left corner 
of the figure. 


to 2 dimensions, (2) uniform surface embedding, and (3) a guided adiabatic 
pull-back of the coordinates from the uniformed surface ptj, to the original 
surface. These three steps ameliorate many difficulties in solving the sys¬ 
tem of equations and provide controllable criteria for obtaining a suitable 
“smooth” solution. We describe these steps in the the next three subections 
and introduce a notion of smoothness. These three subsections of AIM are 
illustrated in hgure. From now on when we mention the metric we are 
referring to a polyhedral metric. 
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2.1 Uniformization via Ricci Flow 


We use a discrete Ricci flow to And a conformal relationship between po and 
Ptf. A conformal factor is assigned to each vertex a of po and is denoted by 
the set {ua{t = 0) = 0}. The relationship between edge lengths of po and pt 
is given by, 


Lb(t) 


Lbi0)e 2 , 


(4) 


where a and b index the bounding vertices of edge iab- 

The discrete Ricci flow equations we use for the conformal factors are 
given by, 

^ (5) 

where ka is the Gaussian curvature at vertex a and k is the target curvature 
of the surface. To keep the surface area constant, we choose k to be the 
average Gaussian curvature over the surface. 


k - 

where x is the Euler characteristic and A is the area of the simplicial surface. 
A detailed description of Gaussian curvature can be seen in appendix 

It was shown by Ghow and Lou that combinatorial Ricci flow and discrete 
Ricci flow are essentially equivalent and that they exponentially converge to a 
surface of constant curvature [^. Therefore given some e <C 1, it is assured 
that for sufficiently long times f/ we can Ricci flow the initial conformal 
factors such that the Gaussian curvature at each vertex is close to fc. 


\\ka{tf)-k\\ 2 <e. (7) 

Once we have {ua{tf)}, we create a continuous family of conformal factors 
between po and p*^ using the linear relation, 

Ua (t) = Ua (tf) . (8) 

We also considered using an exponential relationship between conformal fac¬ 
tors and the set of conformal factors produced at each step during Ricci flow. 
These alternate paths between conformal factors were shown to be computa¬ 
tionally more demanding without any apparent benefit over the linear path. 
We discuss this in more detail in Sec. 12.31 
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2.2 Uniform Surface Embedding 

The uniform surface Embedding is broken into four steps, 


1. Embedding on a sphere 

2. Anneal step A 

3. Convexity Routine 

4. Anneal step B. 


We will now discuss each step in detail and explain their necessity beginning 
with embedding on a sphere. 

Ideally a constant Gaussian curvature surface will lie on a sphere centered 
at the origin with radius, 


r = 


1 

\/k' 


(9) 


This gives an additional constraint, = xf + yf + zf, on the coordinates 
of ptf.We are free to specify the initial line segment iab as {xa,ya,Za} = 


Xb 


Vb 


Zb 


r{r‘^ - Lib) + VLl^r‘^{2r^ - Ll^) 

2r2 

- Lib) - 'jLlbr‘^{‘2r^ - Ll^) 

2r2 

r 

71' 


( 10 ) 

( 11 ) 

( 12 ) 


Given this embedded line segment, we embed vertex c by solving. 


= (Xa - Xcf + {ya - ycf + {Za - Z^f, (13) 

= {.Xb - Xcf + ijjb - ycf + {Zb - Zcf, (14) 

= xl + yl + zf (15) 


There are ordinarily two solutions to these equations which correspond to 
the reflection of vertex c about tab- The translational and rotational degrees 
of freedom are fixed once we embed the three vertices of A^fec on the surface 
of the sphere. 




Figure 2: Triangles made of solid lines have already been embedded and are used as pivots 
to solve for the next set of coordinates. Dotted lines represent edges that are used, in 
conjunction with the radius, to compute the eoordinates of the vertex represented by the 
dots. The blue paint brush lines, or determined edges, are edges that are determined once 
the coordinates of its vertices are found. They are not actually used to compute coordinates. 


Starting with this single embedded triangle Aabc we can embed each of 
the three triangles sharing bounding edges. This procedure is illustrated in 
figure. 1^ In fact given a triangle with a single embedded edge, we embed its 
opposite vertex by solving (13) - (15). This allows us to generate a growing 
network of embedded triangles on the surface of the sphere one triangle at a 
time. The dimensional reduction afforded to us by the Ricci flow essentially 
block-diagonalizes the original sparsely-coupled matrix from the quadratic 
distance equations ([^ into many 3x3 matrices. By looking at shared edges 
between embedded and non-embedded triangles, we can successively embed 
each vertex individually until all vertices a G ptj. lie on the sphere. We choose 
the solution to these quadratic systems such that the dihedral angle between 
edge iab is maximal and less than vr. This solution yields two “unfolded” 
triangles embedded on the surface and is illustrated in hgure. 

Here we explain the necessity and logistics behind anneal step A. The 
embedding-tree procedure will ordinarily lead to inconsistencies with the 
isometric embedding equations ([^. The vertices of a constant curvature 
polyhedron ordinarily will not lie on a sphere of constant radius unless the 
number of vertices Nq is sufficiently large. Inconsistencies occur on the “de¬ 
termined edges”, as illustrated in hgure. whose edge lengths do not agree 
with those in pt^. We correct this disagreement by annealing the coordinates 
on the sphere to the edge lengths of pt^. Annealing is done using Newton’s 
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b 


Figure 3: Similar to figure the dotted lines represent the edges used to compute the 
coordinates of vertex k. The vertices of the solid triangle have already been embedded in 
Lastly, the correct solution is selected by maximizing the distance ik. This is eguivalent 
to maximizing Oab such that it is less than tt . 


method to minimize the L 2 norm of (|^ at time t = tf, 

+ [Va (tf) - Vb {tf)f + [Za (tf) - Zb {tf)]'^Y = 0, 
ab 

( 16 ) 

with the coordinates on the sphere as our initial guess. The surface after 
annealing is uniform in the intrinsic sense but extrinsically it may have neg¬ 
ative curvature in some places. Before we move on to £x this issue, we must 
discuss the isometric variations of polyhedra. 

It is well known that polyhedra are not uniquely determined by their 
edge lengths since this restriction says nothing about their extrinsic curva¬ 
ture. For example, figure shows two polyhedra with identical edge lengths 
and Gaussian curvatures but different extrinsic curvatures. If all the edges 
emanating from a vertex have negative extrinsic curvature, that vertex is 
defined as an inverted vertex. This is seen in figure [^. On average there are 
six edges emanating from a vertex of an arbitrary triangulated polyhedron. 
We hnd it useful to define a “paritally-inverted vertex”. A vertex is defined 
as partially inverted if more than seventy percent of the edges emanating 
from it have negative extrinsic curvature. Given our definitions of inverted 
and partially inverted vertices, we define smoothness. 

Definition. An embedding is considered smooth if it does not contain any 
inverted vertices or partially inverted vertices. 

This dehnition is motivated in part by Nollert and Herlod’s suggestions for 
removing non-smooth depressions from the surface. Through our experi¬ 
ments we found allowing vertices to have thirty percent or less of their edges 
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Figure 4' This is an example of two polyhedron with identical edge lengths and Gaussian 
curvatures but different extrinsic curvatures. 


with negative extrinsic curvature is sufficient for allowing smooth depres¬ 
sions. Not only are these embeddings smooth, they are also the most convex 
embeddings our algorithm can produce. 

An inverted vertex can be “popped out” such that the edges maintain 
their lengths but their extrinsic curvatures are no longer negative. This 
process is local meaning all other vertices maintain their coordinates and all 
other edge lengths remain the same. We show this in figure [^. On the 
other hand a partially inverted vertex can not be removed locally and it 
is not clear what it means for them to be inverted. That is why we need 
the convexity routine during the uniform embedding and GAP procedures 
to apply a sort of outward pressure that makes the embedding as convex 
and smooth as possible. The convexity routine is described as follows: first 
compute the average vector Vavg of the vertices sharing an edge with the 
inverted or partially inverted vertex, second compute the difference vector 
defined as hdi// = v — vZ,g, and third define the new position vector of vertex 
V by Vnew = V + ‘^Vdiff- For an inverted vertex this will get one close to 
the ’’popped out” solution were edge-lengths are invariant. For a partially 
inverted vertex the edge-lengths are no longer the same, but when annealing, 
i.e. minimizing in the space of coordinates, our initial guess is now closer to 
a more convex solution. In both cases annealing after the convexity routine 
returns the edge lengths to their original values within some set tolerance, 
this is the purpose of anneal step B. 
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2.3 Guided Adiabatic Pull-Back (GAP) 


In the final step of the AIM algorithm we pull-back from Pt^ to Pq by finding 
coordinates of Pt for all steps tj G [0,f/]. These coordinates are found by 
minimizing (16) at each tf^j. For tf_i we used the uniform embedding found 
in the previous section as our initial guess. If j = 2, 3, we use the coordinates 
at time as the initial guess for embedding surface p*. For j > 3 we 

extrapolate the coordinates at using the previous three coordinate sets. 
This extrapolation procedure has two purposes: first it brings our initial 
guess closer to the global minimum which decreases the convergence time for 
Newton’s method, second it allows us to take larger steps from tf_j to 
while keeping us within an open ball of the global minimum. This better 
initial guess further decreases the run time of GAP. 

There are several paths that relate conformal factors {ua (tf)} to {ua (0)}. 
We chose the linear path given by (|^. This path is optimal because the 
rate of change of Gaussian curvature is constant at each time step tf^j to 
tf-j-i, which allows one’s step size At = to be constant throughout 

GAP. Here Tsteps is the number of steps taken. If the path was not linear, 
we would need adaptive time stepping to account for the variable change 
in Gaussian curvature. For example, one of the alternative paths that we 
considered was given by the conformal factors produce at each step of Ricci 
flow. Since Ricci flow exponentially uniformizes curvature, there is a greater 
change in curvature near t = 0. Therefore it is necessary to decrease At to 
maintain adiabaticity near this region of rapid change. We also considered 


an exponential path given by the time interval, tj = j where 

a is the parameter that controls the rate of convergence. The exponential 
paths suffers the same problem in adaptive stepping as the Ricci flow path. 
Neither of these paths yielded improved results over the linear path and both 
are computationally more demanding. 

As we mentioned in the previous section, minimization techniques can 
not determine extrinsic curvature. Even after applying annealings and the 
convexity routine to Pt^, it is still possible to begin with an initial surface at 
tf that has many clustered and isolated partial inverted vertices. If nothing 
was done GAP would produce an undesired non-smooth surface. To smooth 
out the embedding we apply the convexity routine at each step. Within 
the space of coordinates, this perturbs our initial guess such that we begin 
near solutions without partially inverted vertices. Although this does not 
immediately remove the problem vertices, it gradually reduces them at each 
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step until they are eliminated. We illustrate this gradual push toward a 
smooth embedding in hgure keeping in mind that (16) is nonlinear thus 
having several global minima whose residuals are zero. Using all of these 


# BAM without convexity routine 

■ BAM with convexity routine 


-► Convexity routine 



Number of inverted and partially inverted vertices 

Figure 5: This is a overly simplified visual representation of minimization in coordinate 
space. Each application of the convexity routine perturbs us away from undesirable global 
minima. After many time steps we eventually reach a global minimum without any inverted 
or partially inverted vertiees. 

elements, we dehned our evolution as adiabatic if At is small enough such 
that we transition to and remain near a global minimum with no inverted or 
partially inverted vertices. This gradual nudging during GAP and its ability 
to sufficiently smooth the surface is the main result of this manuscript. 


3 Numerical Tests 

We present three numerical examples of the AIM algorithm. The hrst two 


example meshes were made using the distmesh program 16 . As input, 


distmesh is given z (a:, y) to produce coordinates and a list of triangles, by 
vertices, for surfaces in We use these coordinates to get edge lengths 
by calculating the Euclidean distance between vertices. We then use the 
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extracted edge lengths together with the list of triangles as the original poly¬ 
hedral metric for which we apply AIM. Once AIM is complete, we check how 
well AIM reproduces intrinsic curvature by comparing edge lengths of the 
distmesh surface to those produced AIM. We compare extrinsic curvature by 
looking at the convergence of integrated mean curvature produced by AIM 
to the continuum value. 

Our third example is a surface just outside the ergosphere of a Kerr 
black hole that is right below maximal rotation. We could not use distmesh 
to triangulate this surface since its embedding equations are quasi-analytic, 
which makes them incompatible with distmesh’s input format. Because of 
this, we triangulated the ergosphere ourselves using the components of the 
metric to assign edge lengths between vertices. 


3.1 Distmesh surfaces 


The hrst two surfaces we chose as our test cases are given by, 

z'^ 

^ + :^ + ^ = {Ellipsoid) 

^ + 71 = 1 (Drum) 


(17) 

(18) 


where a = 3, b = 2, c=l, d = 2, e = l and / = vT~5. Both surfaces 
have strictly positive point wise Gaussian and mean curvatures. We will 
refer to the surfaces described by ( pT| and (18) as the ellipsoid and drum, 
respectively. 

The following work was programmed using matlab. This preliminary ap¬ 
plication of AIM does not attempt to optimize our algorithm. The goal is to 
verify our approach and test how accurately it preserves distances and curva¬ 
ture quantities. Run times can be signihcantly improved by parallelizing the 
code and using preconditioning for our optimization routines. When using 
AIM, one is allowed to choose the accuracy of their embedding by manipulat¬ 
ing the tolerances during the annealing and GAP procedures. Tolerances are 
set using the value of the residuals and the magnitude of the steps in coordi¬ 
nate space for each iteration of the quasi-Newton method. These tolerances 
are set to 10“® and 10“®, respectively. This means the Newton’s method will 
stop if the edge lengths are within 10“® accuracy or if the distance between 
points in coordinate space from one iteration to another is on the order of 
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(a) 


(b) 


(c) 


Figure 6: On the left we have the embedding in using the eontinuum equations. In the 
center we have the smooth embeddings produced by AIM. On the right we have a non¬ 
smooth embeddings which are produced when we do not pull-back adiabatically and we do 
not use the convexity routine. The color coding indicates the mean curvature associated to 
eaeh triangle. The ellipsoid has a resolution of 3226 vertices while the drum’s resolution 
is 3246. 
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10“®. Given these settings, we determine the number of steps necessary for 
maintaining adiabaticity when performing GAP. We also present our analy¬ 
sis on how each portion of the code scales with increased vertices, as well as 
the scaling for the number of time steps necessary for adiabaticity. For each 
surface tested, we recovered an average difference in the edge length between 
the original metric and AIM within our prescribed tolerance. The standard 
deviation in the difference between edge lengths are also on the order of our 
prescribed tolerance. This result was tested with tolerances in accuracy set 
to 10“® and 10“®. Our reported scaling analysis corresponds to a tolerance 
of 10-F 

We test whether or not the extrinsic curvature is recovered by hrst looking 
at a qualitative comparison between continuum and AIM embeddings as seen 
in £gure[^ These surfaces are color coded (shaded) as a function of mean cur¬ 
vature. The mean curvature of the continuum surfaces were computed using 
the continuum equations while the triangulated mean curvatures were com¬ 
puted using discrete methods [^. Figure gives a quantitative comparison, 
and convergence properties, of the mean curvature between the continuum 
and AIM by plotting the percent error in integrated mean curvature as a 
function of resolution. 

We made log-log plots of run time verses number of vertices for each of 
the three sub routines and the AIM algorithm as a whole to understand the 
scaling behavior of our code. The same plots were made for analyzing the 
number of steps necessary for adaibaticity as a function of resolution. The 
plots are given in hgure It is observed that the highest order contribution 
to the scaling of AIM comes from the Ricci flow procedure. Our overall 
results suggest that AIM scales sub cubically. 


3.2 Modified Ergosphere of Kerr spacetime 

In 1973 Larry Smarr analytically embedded an axisymmetric 2-surfaces of a 


rotating black hole geometry 18 . Following this we use the Boyer-Lindquist 


coordinates for the metric representation of the rotating black hole spacetime. 


= gudt^ 2gt^dtd(j) grrdr"^ geedO"^ + 


(19) 
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percent error of integrated mean curvature 



Figure 7 ; Convergence of integrated mean curvature as a function of resolution 


where 


9tt 

9t(t> 

9rr 

9ee 

9(t>(t> 


- 1 - 


2Mr + 


2Mr - 


asin^ 6 , 


S 

A’ 

S, 


2 ^ 2 _|_ (2Mr — sin^ 6 


Here we used the usual definitions where 


sin^ 9. 


S := r^ + a^cos^^, 

A := — 2Mr + + Q^. 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 


(25) 

(26) 


For our third example of the AIM algorithm we embedded a distorted ergo- 
sphere of a nearly maximally rotating black hole geometry with zero charge. 
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Figure 8: These are log-log plots of each section of AIM and the algorithm as a whole. The 
slopes indicate scaling behavior of our algorithm as a function of Nq. Our overall results 
suggest that AIM scales sub cubically with time. Time is in units of seconds. 


This surface is defined on a f = const hypersurface thus making dt = 0. 
We have a/M = 0.975 with a = 0.975 and M = 1 where a/M = 1 gives a 
maximally rotating surface. Given that the surface is axiasymmetric, we can 


write r = R{6) which implies that dr = R^gdO. Inserting this into (19) we 
get the two metric, 


ds^ = S 


Rl 

i + ^]de'^ 


hee 


+ 




sm' 


edej)"^. (27) 
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The ergosphere is defined by the radius in which = 0 in (19). At this 
point the surface is elongated at the poles which makes it an extreme surface 
to embed. To make the poles more smooth we distort the surface by adding 
a small parameter e that gives us a surface slightly outside the ergosphere. 
The radius as a function of 9 is given as, 


r = R (6) = M + ^1— j cos^^j +eM. 


Let the isometric embedding functions be defined as, 

x{9,(j)) = p{9)cos(j), 
y{e,(j)) = p(0)sin0, 

z{e) = f{9). 


Equating the line elements between spaces gives the relation, 
ds^ = dx^ + dy'^ + dz"^ = d9‘^ + p^dcf)^. 

Using the above equations to solve for p and 6 we have. 


p(») = 



(28) 


(29) 

(30) 

(31) 


(32) 


(33) 

(34) 


These equations give us our continuum embedding in The metric in the 
continuum is used to assign edge lengths between vertices for our polyhedral 
metric for which we apply AIM. Figure illustrates the continuum embed¬ 
ding, the smooth embedding given by AIM and the non-smooth embedding. 
The percent error of the integrated mean curvature between the continuum 
and AIM embeddings is 0.62 percent, which comparable to the percent er¬ 
ror for the drum at the same resolution. Although we did not analyze the 
convergence of integrated mean curvature for this surface, we expect its con¬ 
vergence to be similar to the ellipsoid and drum surfaces. Our main obstacle 
in completing this convergence analysis is producing consistent high quality 
triangulations for increasing resolution. The edge lengths are recovered up 
to the set accuracy of 10“®. 
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Figure 9: On the left we have the original embedding made using the continuum embedding 
equations. The middle is the smooth embedding given by AIM. On the right we have the 
non-smooth embedding that is produce when one does not step adiabatically and does not 
use the convexity routine. 
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4 Conclusion 


The AIM algorithm was developed using several techniques from previous 
isometric embedding algorithms. Our tests of AIM included axial and non- 
axial symmetric surfaces with strictly positive Gaussian curvature. We also 
tested a surface just above the ergosphere of a Kerr black hole, which is 
axial-symmetric with regions of negative Gaussian cnrvature. 

The main innovation of the AIM algorithm comes from the guided adia¬ 
batic pull-back. Although we used triangulated surfaces, similar to the wire 
mesh of Nollert and Herold, GAP allows us to distinguish between the global 
minima of smooth and non-smooth polyhedra. GAP also prevents the prob¬ 
lem of getting trapped in local minima, as seen by Bondarescn, Alcnbierre 
and Seidel, during Newton’s method . For the snrfaces tested, we were ca¬ 
pable of reaching residnals on the order of our prescribed tolerance. AIM 
also uses an embedding flow similar to Jasinlek and Korzyhski’s, bnt is not 
restricted to strictly positive scalar curvatnre snrfaces. 

Althongh AIM does not necessarily prodnce the convex polyhedron given 
a convex polyhedral metric, we were capable of dehning and producing a 
smooth embedding. Fnture analysis is planned to provide an upper bound 
of the deviation from convexity for smooth embeddings prodnced by AIM. 
This bonnd should be a function of resolntion and curvature of the surface. 

It was not onr goal or intention to optimize the AIM algorithm. Its rnn 
time and scaling may be signihcantly improved by using preconditioning for 
the optimization rontines and parallelization. We would also like to develop 
an algorithm that prodnces well posed triangnlation, similar to those made 
by distmesh, given any compact two metric. Since our original interest in 
isometric embeddings stemmed from compnting quasi-local energy in general 
relativity, we plan to use AIM for this purpose. 
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A Intrinsic Curvature Calculations 


In Sec. 2.1 we used simplicial Ricci flow to obtain a triangular mesh with con¬ 
stant Gaussian curvature at each vertex. This procedure relied on curvature 
values defined on the surface of our lattice. The purpose of this appendix is 
to inform the reader about the nature of curvature on simplicial geometries 
and provide the exact equation used to compute Gaussian curvature. 

that 


It was established by Alexandrov, and utilized by Tullio Regge 19 


the intrinsic curvature of a simplicial geometry, of arbitrary dimensions, is 
concentrated at co-dimension 2 simplices called hinges. The curvature at 
these hinges is a conic singularity whose value goes to infinity as the area of 
the loop of parallel transport enclosing the hinge goes to zero. Later it was 


determined by W. Miller et. al. 20,21 that given a Daulanay lattice S with 
a circumcentric Voronoi dual lattice 5*, it was possible to define a unique 
path of parallel transport such that the sectional and Gaussian curvatures 
are defined. The area h* of this loop is the Voronoi area perpendicular to a 
hinge in the Daulanay lattice. The Voronoi area h* is constructed using dual 
edges X E S* whose vertices are defined as the circumcenters of d-simplices. 
If one were to parallel transport a vector around this loop, they would find 
their vector rotated by the deficit angle at the hinge. The deficit angle is 
defined as eh = 271 — Vi^ where the sum is over the interior angles of 
d-simplicies that share the hinge h. Using the deficit angle and the dual area 
associated with h, the Gaussian curvature is defined as. 


kh = 


eh 

\h*\ 


(36) 


Everything mentioned in this appendix is applicable for arbitrary dimen¬ 
sions. Since we are embedding an surface, we are only concerned with 
two dimensions. This means the curvature is concentrated on vertices. 
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